stanEx0

options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)

use stan with cmdstanr

install.packages('cmdstanr',repos=c('https://mc-stan.org/r-packages/',getOption('repos')))

library(cmdstanr)
check_cmdstan_toolchain(fix=T)
install_cmdstan(cores=2)
cmdstan_version()


restart R


options(mc.cores=parallel::detectCores())


make stan file '---.stan'

ex0-0.stan

Bernoulli distribution

data{
  int N;
  array[N] int y;
}
parameters{
  real<lower=0,upper=1> p;
}
model{
  p~beta(1,1);
  y~bernoulli(p);
}

execute mcmc sampling function

goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
  mcmc=mdl$sample(
  data=data,
  seed=1,
  chains=4,
  iter_sampling=smp,
  iter_warmup=wrm,
  thin=th,
  refresh=1000
  )
  mcmc
}

data=list(N=9,y=sample(0:1,9,replace=T))

mdl=cmdstan_model('./ex0-0.stan')
#mdl$exe_file()

mle=mdl$optimize(data=data)

## Initial log joint probability = -7.2058 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        5      -6.18265   0.000131667   2.68779e-08           1           1        8    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.3 seconds.

mle

##  variable estimate
##      lp__    -6.18
##      p        0.56

mcmc=mdl$sample(
  data=data,
  #seed=1,
  #chains=4,
  iter_sampling=500,
  iter_warmup=100,
  #thin=1,
  refresh=1000
)

## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.

# or
#mcmc=goMCMC(mdl,data)

mcmc

##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##      lp__ -8.07  -7.80 0.72 0.30 -9.52 -7.58 1.00     1029     1211
##      p     0.55   0.55 0.14 0.14  0.30  0.77 1.00      748     1291

#or
#mcmc$summary()


drw=mcmc$draws('p')

par(mfrow=c(4,1),mar=c(1,5,1,1))
drw[,1,] |> plot(type='l',xlab='',ylab='p')
drw[,2,] |> plot(type='l',xlab='',ylab='p')
drw[,3,] |> plot(type='l',xlab='',ylab='p')
drw[,4,] |> plot(type='l',xlab='',ylab='p')

par(mar=c(3,5,3,3))

par(mfrow=c(1,2))
drw %>% hist()
drw %>% density() %>% plot()

see mcmc result and parameters function

seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
  print(mcmc)
  prs=mcmc$metadata()$model_params[-1] # reject lp__
  for(pr in prs){
    if(grepl('^y',pr)) next # not show predictive value "y*" information
    if(exc!='' && grepl(paste0('^',exc),pr)) next
    drw=mcmc$draws(pr)
    if(ch){
      par(mfrow=c(4,1),mar=c(1,5,1,1))
      drw[,1,] |> plot(type='l',xlab='',ylab=pr)
      drw[,2,] |> plot(type='l',xlab='',ylab=pr)
      drw[,3,] |> plot(type='l',xlab='',ylab=pr)
      drw[,4,] |> plot(type='l',xlab='',ylab=pr)
      par(mar=c(3,5,3,3))
    }

    par(mfrow=c(1,2))
    drw |> hist(main=pr,xlab='')
    drw |> density() |> plot(main=pr)    
  }
}

mcmc$metadata()$stan_variables

## [1] "lp__" "p"

mcmc$metadata()$model_params

## [1] "lp__" "p"

seeMCMC(mcmc)

##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##      lp__ -8.07  -7.80 0.72 0.30 -9.52 -7.58 1.00     1029     1211
##      p     0.55   0.55 0.14 0.14  0.30  0.77 1.00      748     1291

save and load mdl,mcmc result

saveRDS(mdl,'bin_mdl.rds')
saveRDS(mcmc,'bin_mcmc.rds')

mdl=readRDS('bin_mdl.rds')
mcmc=readRDS('bin_mcmc.rds')

library(qs) #faster
qsave(mcmc,'bin_mcmc.qs')
mcmc=qread('bin_mcmc.qs')

use formula and design matrix

n=30
tb=tibble(x=runif(n,0,9),y=rnorm(n,x,1))
tb

## # A tibble: 30 × 2
##        x      y
##    <dbl>  <dbl>
##  1 2.83   2.45 
##  2 4.26   3.76 
##  3 7.60   8.19 
##  4 6.92   7.07 
##  5 5.66   4.49 
##  6 1.60   1.88 
##  7 4.68   3.63 
##  8 8.39  10.3  
##  9 0.615  0.233
## 10 7.94   7.47 
## # ℹ 20 more rows

f0=formula('y~x')
f0

## y ~ x

lm(f0,tb)

## 
## Call:
## lm(formula = f0, data = tb)
## 
## Coefficients:
## (Intercept)            x  
##       -0.29         1.04

X=model.matrix(f0,tb) # get explanatory variables from what in formula
X

##    (Intercept)     x
## 1            1 2.833
## 2            1 4.264
## 3            1 7.599
## 4            1 6.918
## 5            1 5.661
## 6            1 1.600
## 7            1 4.682
## 8            1 8.389
## 9            1 0.615
## 10           1 7.937
## 11           1 7.109
## 12           1 5.520
## 13           1 6.692
## 14           1 8.331
## 15           1 2.101
## 16           1 0.841
## 17           1 0.279
## 18           1 3.446
## 19           1 3.392
## 20           1 0.111
## 21           1 6.372
## 22           1 4.384
## 23           1 8.172
## 24           1 5.698
## 25           1 2.587
## 26           1 3.794
## 27           1 7.949
## 28           1 8.567
## 29           1 7.475
## 30           1 3.166
## attr(,"assign")
## [1] 0 1

tb=tibble(x1=runif(n,0,9),
          x2=runif(n,0,9),
          c1=sample(c('a','b','c'),n,replace=T),
          y=rnorm(n,x1-x2+(c1=='b')*5-(c1=='c')*5,1))
tb

## # A tibble: 30 × 4
##       x1    x2 c1           y
##    <dbl> <dbl> <chr>    <dbl>
##  1 4.34  8.77  b       0.607 
##  2 7.94  7.61  c      -6.19  
##  3 2.28  8.46  c     -10.1   
##  4 3.25  4.62  a       0.0277
##  5 0.743 1.11  b       5.26  
##  6 3.05  6.28  a      -2.81  
##  7 0.371 0.190 b       4.75  
##  8 0.273 0.823 a      -1.39  
##  9 3.65  7.06  a      -3.82  
## 10 5.29  0.805 c      -0.213 
## # ℹ 20 more rows

f0=formula('y~x1')
f0

## y ~ x1

lm(f0,tb)

## 
## Call:
## lm(formula = f0, data = tb)
## 
## Coefficients:
## (Intercept)           x1  
##      -4.809        0.991

model.matrix(f0,tb)

##    (Intercept)    x1
## 1            1 4.341
## 2            1 7.937
## 3            1 2.280
## 4            1 3.253
## 5            1 0.743
## 6            1 3.045
## 7            1 0.371
## 8            1 0.273
## 9            1 3.649
## 10           1 5.294
## 11           1 8.480
## 12           1 1.910
## 13           1 2.685
## 14           1 8.004
## 15           1 5.883
## 16           1 3.892
## 17           1 5.148
## 18           1 8.230
## 19           1 1.437
## 20           1 1.155
## 21           1 5.810
## 22           1 7.283
## 23           1 3.715
## 24           1 0.403
## 25           1 3.792
## 26           1 7.194
## 27           1 8.346
## 28           1 1.072
## 29           1 7.658
## 30           1 8.515
## attr(,"assign")
## [1] 0 1

f1=formula('y~x1+x2')
f1

## y ~ x1 + x2

lm(f1,tb)

## 
## Call:
## lm(formula = f1, data = tb)
## 
## Coefficients:
## (Intercept)           x1           x2  
##       0.909        1.057       -1.380

model.matrix(f1,tb) # get explanatory variables from what in formula

##    (Intercept)    x1    x2
## 1            1 4.341 8.775
## 2            1 7.937 7.608
## 3            1 2.280 8.455
## 4            1 3.253 4.624
## 5            1 0.743 1.110
## 6            1 3.045 6.284
## 7            1 0.371 0.190
## 8            1 0.273 0.823
## 9            1 3.649 7.058
## 10           1 5.294 0.805
## 11           1 8.480 6.943
## 12           1 1.910 1.387
## 13           1 2.685 8.230
## 14           1 8.004 2.346
## 15           1 5.883 2.871
## 16           1 3.892 6.618
## 17           1 5.148 5.041
## 18           1 8.230 1.213
## 19           1 1.437 3.285
## 20           1 1.155 5.542
## 21           1 5.810 4.366
## 22           1 7.283 2.415
## 23           1 3.715 4.136
## 24           1 0.403 1.211
## 25           1 3.792 5.762
## 26           1 7.194 5.203
## 27           1 8.346 1.504
## 28           1 1.072 7.879
## 29           1 7.658 6.016
## 30           1 8.515 2.861
## attr(,"assign")
## [1] 0 1 2

f2=formula('y~x1+x2+c1')
f2

## y ~ x1 + x2 + c1

lm(f2,tb)

## 
## Call:
## lm(formula = f2, data = tb)
## 
## Coefficients:
## (Intercept)           x1           x2          c1b          c1c  
##      -0.102        1.088       -0.969        5.046       -5.644

model.matrix(f2,tb) # make categorical variable to dummy variable

##    (Intercept)    x1    x2 c1b c1c
## 1            1 4.341 8.775   1   0
## 2            1 7.937 7.608   0   1
## 3            1 2.280 8.455   0   1
## 4            1 3.253 4.624   0   0
## 5            1 0.743 1.110   1   0
## 6            1 3.045 6.284   0   0
## 7            1 0.371 0.190   1   0
## 8            1 0.273 0.823   0   0
## 9            1 3.649 7.058   0   0
## 10           1 5.294 0.805   0   1
## 11           1 8.480 6.943   0   1
## 12           1 1.910 1.387   1   0
## 13           1 2.685 8.230   0   0
## 14           1 8.004 2.346   0   0
## 15           1 5.883 2.871   0   1
## 16           1 3.892 6.618   0   1
## 17           1 5.148 5.041   1   0
## 18           1 8.230 1.213   0   1
## 19           1 1.437 3.285   0   1
## 20           1 1.155 5.542   0   1
## 21           1 5.810 4.366   1   0
## 22           1 7.283 2.415   1   0
## 23           1 3.715 4.136   0   0
## 24           1 0.403 1.211   0   1
## 25           1 3.792 5.762   0   0
## 26           1 7.194 5.203   0   0
## 27           1 8.346 1.504   1   0
## 28           1 1.072 7.879   0   1
## 29           1 7.658 6.016   0   1
## 30           1 8.515 2.861   0   0
## attr(,"assign")
## [1] 0 1 2 3 3
## attr(,"contrasts")
## attr(,"contrasts")$c1
## [1] "contr.treatment"

use brms

library(brms)

rst=brm(f2,
        data=tb,
        family=gaussian('identity'), # binomial('logit'),poisson('log')
        #seed=1,
        #chains=4,
        iter=500,
        warmup=100,
        #thin=1,
        prior=c(set_prior('',class='Intercept'),
                set_prior('',class='sigma'))
        )

rst

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: y ~ x1 + x2 + c1 
##    Data: tb (Number of observations: 30) 
##   Draws: 4 chains, each with iter = 500; warmup = 100; thin = 1;
##          total post-warmup draws = 1600
## 
## Regression Coefficients:
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.08      0.48    -1.04     0.86 1.00      894     1070
## x1            1.09      0.06     0.98     1.20 1.00     1946     1096
## x2           -0.97      0.06    -1.10    -0.85 1.00     1596     1348
## c1b           5.02      0.44     4.19     5.90 1.01      545      708
## c1c          -5.64      0.38    -6.40    -4.93 1.00      605      612
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.86      0.12     0.66     1.13 1.00     1172     1134
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

stancode(rst)

## // generated with brms 2.21.0
## functions {
## }
## data {
##   int<lower=1> N;  // total number of observations
##   vector[N] Y;  // response variable
##   int<lower=1> K;  // number of population-level effects
##   matrix[N, K] X;  // population-level design matrix
##   int<lower=1> Kc;  // number of population-level effects after centering
##   int prior_only;  // should the likelihood be ignored?
## }
## transformed data {
##   matrix[N, Kc] Xc;  // centered version of X without an intercept
##   vector[Kc] means_X;  // column means of X before centering
##   for (i in 2:K) {
##     means_X[i - 1] = mean(X[, i]);
##     Xc[, i - 1] = X[, i] - means_X[i - 1];
##   }
## }
## parameters {
##   vector[Kc] b;  // regression coefficients
##   real Intercept;  // temporary intercept for centered predictors
##   real<lower=0> sigma;  // dispersion parameter
## }
## transformed parameters {
##   real lprior = 0;  // prior contributions to the log posterior
## }
## model {
##   // likelihood including constants
##   if (!prior_only) {
##     target += normal_id_glm_lpdf(Y | Xc, Intercept, b, sigma);
##   }
##   // priors including constants
##   target += lprior;
## }
## generated quantities {
##   // actual population-level intercept
##   real b_Intercept = Intercept - dot_product(means_X, b);
## }

standata(rst)

## $N
## [1] 30
## 
## $Y
##  [1]   0.6071  -6.1854 -10.1209   0.0277   5.2596  -2.8055   4.7511  -1.3912
##  [9]  -3.8195  -0.2133  -2.0143   5.8267  -5.3684   6.5560  -1.4250  -8.8055
## [17]   6.1751   1.6742  -7.2418 -10.0584   8.1638  10.1839  -0.2569  -6.6813
## [25]  -0.6644   3.0950  11.4923 -13.1887  -2.9032   5.7330
## 
## $K
## [1] 5
## 
## $Kc
## [1] 4
## 
## $X
##    Intercept    x1    x2 c1b c1c
## 1          1 4.341 8.775   1   0
## 2          1 7.937 7.608   0   1
## 3          1 2.280 8.455   0   1
## 4          1 3.253 4.624   0   0
## 5          1 0.743 1.110   1   0
## 6          1 3.045 6.284   0   0
## 7          1 0.371 0.190   1   0
## 8          1 0.273 0.823   0   0
## 9          1 3.649 7.058   0   0
## 10         1 5.294 0.805   0   1
## 11         1 8.480 6.943   0   1
## 12         1 1.910 1.387   1   0
## 13         1 2.685 8.230   0   0
## 14         1 8.004 2.346   0   0
## 15         1 5.883 2.871   0   1
## 16         1 3.892 6.618   0   1
## 17         1 5.148 5.041   1   0
## 18         1 8.230 1.213   0   1
## 19         1 1.437 3.285   0   1
## 20         1 1.155 5.542   0   1
## 21         1 5.810 4.366   1   0
## 22         1 7.283 2.415   1   0
## 23         1 3.715 4.136   0   0
## 24         1 0.403 1.211   0   1
## 25         1 3.792 5.762   0   0
## 26         1 7.194 5.203   0   0
## 27         1 8.346 1.504   1   0
## 28         1 1.072 7.879   0   1
## 29         1 7.658 6.016   0   1
## 30         1 8.515 2.861   0   0
## attr(,"assign")
## [1] 0 1 2 3 3
## attr(,"contrasts")
## attr(,"contrasts")$c1
## [1] "contr.treatment"
## 
## 
## $prior_only
## [1] 0
## 
## attr(,"class")
## [1] "standata" "list"

plot(rst) #mcmc trace plot

mcmc=as.mcmc(rst,combine_chains=T) #mcmc sampling

x_new=tibble(x1=runif(3,0,9),x2=runif(3,0,9),c1=c('a','b','c'))
x_new

## # A tibble: 3 × 3
##      x1    x2 c1   
##   <dbl> <dbl> <chr>
## 1 3.98   2.88 a    
## 2 0.862  2.32 b    
## 3 1.13   8.01 c

# bayes credible interval
fitted(rst,x_new)

##      Estimate Est.Error    Q2.5  Q97.5
## [1,]     1.45     0.309   0.844   2.07
## [2,]     3.63     0.368   2.910   4.37
## [3,]   -12.28     0.371 -13.010 -11.55

marginal_effects(rst,effects='x1') |> plot(points=T)

marginal_effects(rst,effects='x2') |> plot(points=T)

marginal_effects(rst,effects='x1:c1') |> plot(points=T)

marginal_effects(rst,effects='x2:c1') |> plot(points=T)

# bayes predicted interval
predict(rst,x_new)

##      Estimate Est.Error    Q2.5  Q97.5
## [1,]     1.47     0.927  -0.338   3.27
## [2,]     3.60     0.953   1.641   5.57
## [3,]   -12.27     0.944 -14.130 -10.32

marginal_effects(rst,effects='x1',method='predict') |> plot(points=T)

marginal_effects(rst,effects='x2',method='predict') |> plot(points=T)

marginal_effects(rst,effects='x1:c1',method='predict') |> plot(points=T)

marginal_effects(rst,effects='x2:c1',method='predict') |> plot(points=T)

stan syntax

code block
  add semicolon ";" at the end of sentence

  data{
    define objects from R as list
  }
  parameters{
    define estimating parameters from MCMC, can't use int type
  }
  transformed parameters{
    define caluculating parameters from other parameters
  }
  model{
    caluculate log liklihood
    define stochastic model, ex. y~normal(m,s) or target+=dist._lpdf(y|m,s)
  }
  generated quantities{
    define caluculating quantities from parameters
  }


object type and its difinition of single and array

  variable
    integer int x, array[N] int x, array[N,M] int x
    real    real x, array[N] real x, array[N,M] real x
      for array[N] int/real x
        x[i]      x_i
        x[i1:i2]  array (x_i1...x_i2)
      for array[N,M] int/real x
        x[i,j]    x_ij
        x[i]      array (x_i1...x_iM)  
  
  vector
    vector[K] x, array[N] vector[K] x
      for vector[K] x
        x[i]      x_i
        x[i1:i2]  vector (x_i1...x_i2)
      for array[N] vector[K] x
        x[i,k]    x_ik
        x[i]      vector[k] x[i]
  
    row_vector[K] x, array[N] row_vector[K] x
    simplex[K] x        x_i[0,1], Σ x_i=1
    unit_vector[K] x    Σ x_i^2=1
    ordered[K] x        x_1 < x_2 ... < x_K
  
  matrix
    matrix[J,K] x, array[N] matrix[J,k] x
    cov_matrix[K] x     symmetric, all eigen>=0
    corr_matrix[K] x    symmetric, all eigen>=0, x_ij [0,1], x_ii=1
      for matrix[J,K] x
        x[j,k]    x_jk
        x[j,]     row_vector[K] x_j
        x[,k]     vector[J] x_k
        x[j1:j2,k1:k2]  matrix[j2-j1+1,k2-k1+1]
      for array[N] matrix[J,K] x
        x[i]      matrix[J,K] x_i

    to fasten a loop for doing matrix[J,K] x
        for(k in 1:K){
          for(j in 1:J){
            do x[j,k]
          }
        }


parameter constraint

  real<lower=0,upper=1> x;
  real<lower=l,upper=u> x;
  
  real a[N];
  real<lower=min(a),upper=max(a)> x;
  
  vector<lower=0,upper=1>[K] x;

using simplex
  categorical distribution
    y~Cat(h)  
      i=1~n, k=1~K, y_i[1,K], h_k[0,1], Σ h_k=1
      P[y=k]=h_k

  data{
    int N;
    int K;
    array[N] int<lower=1,upper=K> y; // R y[n],y_i[1,K]
  }
  parameters{
    simplex[K] h;
  }
  model{
    for(i in 1:N){
      y[i]~categorical(h);
    }
  }

  multinomial distribution
    y~multi(n,h)
      i=1~n, k=1~K, c_i[1:K], y_k=#(c_i=k), y_k[0,n]  
      h_k[0,1], Σ h_k=1  
      P[c=k]=h_k  

  data{
    int K;
    array[K] int<lower=0> y; // R y=table(factor(c,levels=1:K))
  }
  parameters{
    simplex[K] h;
  }
  model{
    y~multinomial(h);
  }


use vector, matrix instead of array to be fast

  using array
    data{
      int N;
      array[N] real x;
      array[N] real y;
    }
    parameters{
      real<lower=0> s;
    }
    transformed parameters{
      array[N] real m;
      for(i in 1:N){
        m[i]= ; // caluculation of x[i]
      }
    }
    model{
      for(i in 1:N){
        y[i]~normal(m[i],s);
    }

  using vector
    data{
      int N;
      vector[N] x;
      vector[N] y;
    }
    parameters{
      real<lower=0> s;
    }
    transformed parameters{
      vector[N] m;
      m= ; //caluculation of x
    }
    model{
      y~normal(m,s);
        or
      target+=normal_lpdf(y|m,s)
    }


    mixed model
      data{
        int N;
        int K;
        array[N] int<lower=1,upper=K> ID;
        vector[N] x;
        vector[N] y;
      }
      parameters{
        real a0;
        real b0;
        vector[K] a;
        vector[K] b;
        real<lower=0> sa;
        real<lower=0> sb;
        real<lower=0> s;
      }
      model{
        a~normal(a0,sa);
        b~normal(b0,sb);
        y~normal(a[ID]+b[ID].*x,s); // .* is multiply element by element
      }
    
    
    K dimension multi normal distribution
      data{
        int N;
        int K;
        arrat[N] vector[K] y;
      }
      parameters{
        vector[K] m;
        cov_matrix[K] cov;
      }
      model{
        y~multi_normal(m,cov);
      }
    
    
    multi regression
      data{
        int N;
        int K;
        matrix[N,K] x; // model.matrix x from R 
        vector[N] y;
      }
      parameters{
        vector[k] b;
        real<lower=0> s;
      }
      transformed parameters{
        vector[k] m;
        m=x*b;
      }
      model{
        y~normal(m,s);
      }

ex0-1.stan

categorical distribution from raw value

data{
  int N;
  int K;
  array[N] int<lower=1,upper=K> y; // R y[n],yi[1:K]
}  
parameters{
  simplex[K] h;
}
model{
  y~categorical(h);
}

c0=c(1,2,3,4)
h=c(0.4,0.3,0.2,0.1)
y=sample(c0,20,h,replace=T)
table(y) |> prop.table()

## y
##    1    2    3    4 
## 0.35 0.15 0.45 0.05

data=list(N=length(y),K=length(c0),y=y)

mdl=cmdstan_model('./ex0-1.stan')

mle=mdl$optimize(data=data)

## Initial log joint probability = -37.0527 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       10      -23.2224   0.000411501   1.05642e-05           1           1       13    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.

mle

##  variable estimate
##      lp__   -23.22
##      h[1]     0.35
##      h[2]     0.15
##      h[3]     0.45
##      h[4]     0.05

system.time({
  mcmc=goMCMC(mdl,data)
  seeMCMC(mcmc,ch=T)
})

## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.4 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -31.28 -30.89 1.34 1.03 -34.09 -29.87 1.00      891      972
##      h[1]   0.33   0.33 0.09 0.09   0.19   0.48 1.00     3305     1692
##      h[2]   0.17   0.16 0.08 0.08   0.06   0.31 1.00     1524     1090
##      h[3]   0.42   0.42 0.10 0.10   0.26   0.58 1.00     1376     1474
##      h[4]   0.08   0.07 0.06 0.05   0.01   0.20 1.00      819      659

##   ユーザ システム     経過 
##    1.003    0.259    1.290

ex0-2.stan

categorical distribution from frequency

data{
  int K;
  array[K] int<lower=0> y; // R y=table(factor(c,levels=1:K))
}
parameters{
  simplex[K] h;
}
model{
  y~multinomial(h);
}

c0=c(1,2,3,4)
h=c(0.4,0.3,0.2,0.1)
c=sample(c0,20,h,replace=T)
table(c) |> prop.table()

## c
##    1    2    3    4 
## 0.35 0.30 0.15 0.20

y=table(factor(c,levels=1:length(c0)))

data=list(K=length(c0),y=y)

mdl=cmdstan_model('./ex0-2.stan')

mle=mdl$optimize(data=data)

## Initial log joint probability = -30.3904 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        8      -26.7017   0.000171529   0.000109974      0.9636      0.9636       11    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.

mle

##  variable estimate
##      lp__   -26.70
##      h[1]     0.35
##      h[2]     0.30
##      h[3]     0.15
##      h[4]     0.20

system.time({
  mcmc=goMCMC(mdl,data)
  seeMCMC(mcmc,ch=T)
})

## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -33.98 -33.66 1.25 1.06 -36.36 -32.61 1.00     1142     1543
##      h[1]   0.34   0.33 0.10 0.10   0.18   0.50 1.00     3089     1743
##      h[2]   0.29   0.28 0.09 0.09   0.15   0.45 1.01     2464     1680
##      h[3]   0.17   0.16 0.07 0.07   0.06   0.30 1.00     1258     1422
##      h[4]   0.21   0.20 0.08 0.08   0.08   0.36 1.00     1575     1682

##   ユーザ システム     経過 
##    0.907    0.229    1.192

WAIC

ex0-3.stan

data{
  int N;
  vector[N] x;
  vector[N] y;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  target+=normal_lpdf(y|b0+b1*x,s);
}
generated quantities{
  vector[N] ll;
  for(i in 1:N){
    ll[i]=normal_lpdf(y[i]|b0+b1*x[i],s);
  }
}

n=10
x=runif(n,0,9)
y=rnorm(n,x,1)
data=list(N=n,x=x,y=y)

mdl=cmdstan_model('./ex0-3.stan')

mle=mdl$optimize(data=data)

## Initial log joint probability = -30.4205 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       17       -13.284   0.000308314   0.000853684      0.9301      0.9301       24    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.

mle

##  variable estimate
##    lp__     -13.28
##    b0        -0.02
##    b1         0.98
##    s          0.91
##    ll[1]     -1.84
##    ll[2]     -1.04
##    ll[3]     -1.20
##    ll[4]     -0.90
##    ll[5]     -1.06
##    ll[6]     -1.16
##    ll[7]     -0.85
##    ll[8]     -1.53
##    ll[9]     -1.11
##    ll[10]    -2.59

system.time({
  mcmc=goMCMC(mdl,data)
  seeMCMC(mcmc,exc='ll',ch=F)
})

## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -15.16 -14.82 1.45 1.24 -17.90 -13.55 1.01      479      536
##    b0      -0.07  -0.03 0.79 0.68  -1.35   1.09 1.01      384      312
##    b1       0.99   0.98 0.15 0.14   0.76   1.23 1.00      426      390
##    s        1.23   1.15 0.39 0.32   0.76   1.98 1.00      752      756
##    ll[1]   -1.93  -1.79 0.64 0.51  -3.22  -1.16 1.00     1265     1409
##    ll[2]   -1.34  -1.31 0.35 0.34  -1.97  -0.82 1.00      599      635
##    ll[3]   -1.47  -1.41 0.43 0.37  -2.23  -0.91 1.00      954     1220
##    ll[4]   -1.22  -1.20 0.31 0.31  -1.77  -0.76 1.00      615      633
##    ll[5]   -1.30  -1.28 0.30 0.29  -1.83  -0.85 1.00     1010      982
##    ll[6]   -1.39  -1.38 0.32 0.31  -1.98  -0.91 1.00      701      865
##    ll[7]   -1.16  -1.14 0.30 0.31  -1.69  -0.72 1.00      716      760
##    ll[8]   -1.74  -1.61 0.63 0.52  -2.95  -0.98 1.01      470      395
##    ll[9]   -1.34  -1.31 0.30 0.30  -1.90  -0.90 1.00      679      676
##    ll[10]  -2.42  -2.22 0.86 0.72  -4.05  -1.37 1.00     1162     1248

##   ユーザ システム     経過 
##    0.579    0.236    0.831

ll=mcmc$draws('ll') |>
  posterior::as_draws_df() |> 
  select(contains('ll'))
lppd=sum(log(colMeans(exp(ll))))
pwaic=sum(apply(ll,2,var))
waic=-2*lppd+2*pwaic

ex0-4.stan

data {
  int N;
  vector[N] x;
  array[N] int y;
}
parameters {
  real b0;
  real b1;
}
model {
  target+=poisson_lpmf(y | exp(b0+b1*x));
}
generated quantities {
  vector[N] ll;
  for (i in 1:N) {
    ll[i] = poisson_lpmf(y[i] | exp(b0+b1*x[i]));
  }
}

n=10
x=runif(n,-1,1)
y=rpois(n,exp(x))
data=list(N=n,x=x,y=y)

mdl=cmdstan_model('./ex0-4.stan')

mle=mdl$optimize(data=data)

## Initial log joint probability = -52.252 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       10      -10.0548    0.00250835   0.000147357           1           1       12    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.

mle

##  variable estimate
##    lp__     -10.05
##    b0        -0.30
##    b1         1.25
##    ll[1]     -1.22
##    ll[2]     -0.39
##    ll[3]     -0.42
##    ll[4]     -1.35
##    ll[5]     -1.13
##    ll[6]     -1.20
##    ll[7]     -0.25
##    ll[8]     -0.98
##    ll[9]     -1.49
##    ll[10]    -1.61

system.time({
  mcmc=goMCMC(mdl,data)
  seeMCMC(mcmc,exc='ll',ch=F)
})

## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -11.04 -10.76 0.99 0.71 -12.92 -10.10 1.00      820      802
##    b0      -0.44  -0.40 0.42 0.41  -1.16   0.22 1.00      984     1048
##    b1       1.31   1.29 0.61 0.60   0.38   2.29 1.01      692      620
##    ll[1]   -1.38  -1.29 0.34 0.30  -2.06  -1.01 1.00      757      871
##    ll[2]   -0.40  -0.35 0.24 0.21  -0.87  -0.11 1.01      734      798
##    ll[3]   -0.42  -0.37 0.25 0.22  -0.90  -0.12 1.01      744      798
##    ll[4]   -1.50  -1.40 0.25 0.13  -2.04  -1.31 1.00     1557     1343
##    ll[5]   -1.26  -1.18 0.26 0.21  -1.78  -1.00 1.00      796      915
##    ll[6]   -1.27  -1.15 0.31 0.20  -1.94  -1.00 1.00     1749     1450
##    ll[7]   -0.29  -0.22 0.23 0.17  -0.71  -0.05 1.01      697      661
##    ll[8]   -0.93  -0.89 0.33 0.31  -1.56  -0.47 1.00     1375     1412
##    ll[9]   -1.70  -1.59 0.55 0.52  -2.75  -1.03 1.01      710      724
##    ll[10]  -1.90  -1.71 0.53 0.29  -3.01  -1.50 1.00     1242     1295

##   ユーザ システム     経過 
##    0.460    0.226    0.734

ll=mcmc$draws('ll') |>
  posterior::as_draws_df() |> 
  select(contains('ll'))
lppd=sum(log(colMeans(exp(ll))))
pwaic=sum(apply(ll,2,var))
waic=-2*lppd+2*pwaic